Screening of Coulomb Impurities in Graphene 
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We calculate exactly the vacuum polarization charge density in the field of a subcritical Coulomb 
impurity, Z\e\/r, in graphene. Our analysis is based on the exact electron Green's function, obtained 
by using the operator method, and leads to results that are exact in the parameter Za, where a is 
the "fine structure constant" of graphene. Taking into account also electron-electron interactions in 
the Hartree approximation, we solve the problem self-consistently in the subcritical regime, where 
the impurity has an effective charge -Zgfj, determined by the localized induced charge. We find that 
an impurity with bare charge Z = 1 remains subcritical, Z e ^a < 1/2, for any a, while impurities 
with Z = 2, 3 and higher can become supercritical at certain values of a. 
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It has been known for a long time that the single elec- 
tron dynamics in a monolayer of graphite (graphene) is 
described by a massless two-component Dirac equation 
[3, 0- S] ■ A surge of interest in the problem was caused 
by the recent successful fabrication of graphene [4j and 
measurements of transport properties 0, @, 0, HI, M fioj ] , 
including an unconventional form of the quantum Hall ef- 
fect. Due to the Coulomb interaction between electrons, 
graphene represents a peculiar two-dimensional (2D)ver- 
sion of massless Quantum Electrodynamics (QED) [3[. It 
appears to be much simpler than conventional QED be- 
cause the interaction is described by the instantaneous 
1/r Coulomb's law. On the other hand the Fermi veloc- 
ity vf ~ 10 6 m/s ~ c/300 (c is the velocity of light), and 



therefore the "fine structure constant" a 



2 /Hv F 



leading to a strong-coupling version of QED. Below we 
set h = vp = 1. Screening of a charged nucleus due 
to vacuum polarization is an effect of fundamental im- 
portance in QED. This problem was investigated in de- 
tail both in the subcritical and supercritical regimes 
11 , 12, HI, | ■ The problem of charged impurity screen- 
ing in graphene, which also can be treated in terms of 
vacuum polarization, has recently received a lot of atten- 
tion [H, [3, E3, 03, Gj, H3, HE S Hi , due to the impor- 
tance of the problem for transport properties involving 
charged impurities, as well as for our general understand- 
ing of the theory of graphene. 

To leading order in the weak coupling expansion, 
Za <C 1, the induced charge is negative and localized at 
the impurity position, pi nc i = — \e\^(Za)5(r) : which leads 
to screening of the impurity potential 20, 21, 22, 24 1. We 
denote by Z\e\ the impurity charge, and e = — |e| is the 
effective electron charge; from now on we refer to Z as the 
impurity charge with the understanding that it is mea- 
sured in units of |e|. In graphene, the strong-coupling 
problem Za ~ 1 was recently addressed (20j . and it was 
found that the supercritical regime occurs for Za > 1/2, 
where a 1/r 2 tail appears in the induced charge density, 
while in the subcritical regime Za < 1/2, the induced 



charge is always localized at the impurity site. Analyt- 
ical results were also supplemented by numerical lattice 
calculations 2l|, leading to similar conclusions. The in- 
duced charge behavior for small Za was first emphasized 
in a different context in Ref. (24| . and it is in agreement 
with a recent pcrturbative calculation of non-linear vac- 
uum polarization at order (Za) 3 It was also pointed 
out [2 21 ] that a power law tail can appear even in the sub- 
critical regime due to interaction effects (at order Za 2 ). 

In the present paper we investigate analytically the 
induced charge density in graphene in the subcritical 
regime. We use the method of calculation suggested in 
Ref. [IH , where the the induced charge density in a strong 
Coulomb field was obtained in coordinate space in three- 
dimensional (3D) QED. We express the induced charge 
density via the exact Green's function in a Coulomb field, 
calculated within the operator technique [25| . Our main 
result is an exact expression for the polarization charge, 
non-perturbative in the parameter Za, and we explore 
its physical consequences. The exact result allows us to 
determine the effective impurity charge Z & q in a self- 
consistent way. Perhaps most surprisingly, we find that 
an impurity with bare charge Z — 1 can never become 
supercritical, i.e. Z Q ga < 1/2. In spite of this, screening 
can be very substantial, i.e. one can have <C 1 for 
large enough a. 

Electron Green's function in a Coulomb field. It is con- 
venient, for technical reasons, to introduce a small "elec- 
tron mass" M which we set to zero at the end of our 
calculations. This will allow us to avoid some difficulties 
which appear in the calculation of the induced charge 
in massless QED [12( , and will serve the renormalization 
purpose. Physically the mass describes a small energy 
splitting between carbon atoms in the unit cell. We will 
consider only the half-filled case of graphene, i.e. we set 
the chemical potential to zero. 

The electron Green's function G(r, r'|e) in a Coulomb 
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field satisfies the two-component equation 

+■ ^ - (a ■ p) - <t 3 m) G(r, r'|e) = 5(r - r'). (1) 

Here a = (<7 l5 er 2 ), and 01,2,3 are the Pauli matrices; p = 
(Px,P y ) is the momentum operator. Following Ref. [25 1 
we represent the solution of Eq. (fTJ) in the form 



G(r,r'|e) 



ds 
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exp ( is[rA r + — + r(e 2 - M 2 )] 



xJ-6(r-r')5(<l>-<f>'), 



(2) 



where A r is the radial part of the Laplace operator and 
K = d 2 /d(/) 2 + (Za) 2 - iZa(a- n), n = r/r. Then we 
express 8(<fr — </>') via the projectors Pa that are eigen- 
functions of the operator K, KP\ = —\ 2 P\: 



(3) 



x = 7 T 2 . 7 = V* 2 - {Za) 2 , x = m + 1/2 , 

and m = 0, 1, 2, .... The explicit form of Pa for A = 7 — | 
is 



1 /Pi Pi: 



Pa( ^ ,} = 4^ iKi P22 
P11 = (7 + K)e im ^-*"> + (7 - ^) e -' i ( m+1 )^- < >'), 
P12 = -iZa (V 4 ^ 1 )^™^' + e ™* e " i(m+iw ') , 

Pi = ^22! ^*12 = — ^21- 

The projector with eigenvalue A = 7 + 1/2 is ob- 
tained from ((4]) by replacing 7 — > —7. After sub- 
stitution of the expansion ([3]) into Eq. ([2]), the prob- 
lem is reduced to calculation of the action of the op- 
erator exp (— 2is [Ai + & 2 4r]) on ^pS(r — r'). Here 
k 2 = M 2 — e 2 and the operators are defined as: A\ = 
5 (~BF r TFr + ^),A 3 = r. Along with A 2 = -i(rd/dr+ 
1/2), these operators generate an 0(2, 1) algebra, which 
was considered in Ref. [HI in relation to the 3D Green's 
function for the Dirac equation of an electron in a 
Coulomb field. Therefore we can directly use the op- 
erator transformation of that work. As a result we find 
the following integral representation for the solution of 
Eq. Q: 

G(r,v'\e) = -(e+^ + (a-p) + a 3 M^ (5) 



x exp [ik(r + r') cot(fcs) — iirX] J 2 \ 



2fc\/rr' 
sin(fcs) 



Here J„(x) is the Besscl function. 

Induced charge. The induced charge density is 



Pmd(r) = -ieN 



'c 



^Tr{G(r,r|e 



-m, (6) 



where G is the Green's function calculated above, and 
N = 4 reflects the spin and valley degeneracies. The 
contour of integration C goes below the real axis in the 
left half-plane and above the real axis in the right half- 
plane of the complex e plane. Taking into account the 
analytical properties of the Green's function, the contour 
of integration with respect to e can be deformed to co- 
incide with the imaginary axis. The integration contour 
with respect to s in Eq. ([5]) can then also be rotated 
to coincide with the imaginary axis so that it extends 
from zero to — ioo for Ime > 0, and from zero to zoo for 
Ime < 0. After these transformations and obvious change 
of variables, we obtain: 



Pmd{r) = I J deds 



D — y cosh s 



(7) 



x \ 2Zacos(p,s) coths/ 27 (j/) - sin(^s)-y/ 27 (y) 



(4) where k = \/? 



M 2 



2rfc/sinhs, fi = 2Zae/k, 



Ij(x) is the modified Bcssel function of the first kind, 
and 7 27 = dl 2l (y)/dy. We note that p m d{r)-, Eq. 0, is 
an odd function of Za. 

In fact, the expression |J7J) is not well defined and the 
answer depends on the order of integration over e and s. 
To overcome this problem we follow the usual procedure 
of QED. We perform the regularization of the integrals 
introducing a finite upper limit of integration for e, and 
a lower limit of integration for s. Then we carry out 
the renormalization using the obvious physical require- 
ment of zero total induced charge. We can satisfy this 
requirement due to the non-zero mass M because the 
induced charge density diminishes rapidly at distances 
r S> 1/M. After renormalization the result is indepen- 
dent of the cutoff parameters and the order of integration. 
It is technically convenient to do the renormalization in 
momentum space. All details of calculations which we 
need to perform are similar to those described in de- 
tail in Ref. 13[ for the problem of 3D vacuum polariza- 



tion in QED. Finally, we obtain the renormalizcd induced 
charge density in momentum representation pfL d (<i/M). 
The leading term of the asymptotics of this function at 
M — > (or q/M — > 00) is a constant, Qmd- Therefore 
the induced charge density in coordinate space has the 
form 



ici<5(r) + Pdistr 



(8) 
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FIG. 1: The induced charge Qi„d/e = A as a function of Za. 
The solid line is the exact result from Eq. (|9]) and the dashed 
line is the one loop result ,4 (1) = 7r(Za)/2. 



For the induced charge Qind we find 

Qind = eN \^Za + A(Za)] = -\e\A(Za) 
2 °° 

A(Za) = - Im 



m=0 



1 



lnr(7 — iZa) + — ln(7 — iZa) 



iZa 



(7 — iZa)ip(^ — iZa) H iZaxip' (x) 

2x 



(9) 



where is the gamma function and ip(x) = 

d\nT(x)/dx. The induced charge Qind is negative (the 
function A(Za) is positive, see below.) The distributed 
charge density pdistr in Eq. ([8|) is positive and at distances 
r < l/M, p distr (x M 2 ln(l/Mr). The density van- 
ishes in the limit M — > 0, however the total distributed 
charge is —Qind- We comment that in conventional 3D 
QED, the induced charge density also consists of local 
and distributed parts. In that case the distributed den- 
sity Pdistr oc — \e\Za/r 3 at r <C l/M [26| . Interestingly, 
the signs of the local and distributed charges in 3D QED 
arc different from those in Eq. ([8]). The difference is in 
the leading (linear in Za) contribution, while the signs 
of the next-to-leading contributions (non-linear vacuum 
polarization) are the same in both cases. Physically this 
difference is related to the fact that the effective charge e 2 
is renormalizcd (increases at short distances) in conven- 
tional QED, while in graphene the charge is not renor- 
malized 

The series expansion of the function A(Za) = 
(n/2)Za + 4A(Za), for small Za reads 



A(Za) = -(Za) 



0.783(Za) 3 + 1.398(Za) E 



(10) 



The first term of the expansion in Eq. (flTjl) reflects the 
linear one-loop polarization contribution [2 (J . \22\ . Our 
coefficient of the (Za) 3 term has an opposite sign with 
respect to the one found recently in Ref. [22ft . 




FIG. 2: (Color online.) The effective coupling Z e ^a from 
Eq. (fT2)l as a function of Za for Z = 1 (black line), Z = 2 
(red line), and Z — 3 (blue line). Inset: Reduction of the 
impurity coupling relative to its initial value, Z e g/Z. 



The coefficients in the non-linear terms are not small 
and grow as the order increases; this is in contrast with 
conventional QED where they are very small and de- 
crease rapidly 12j. The exact function A(Za) is shown 



in Fig. Q] For Za > 0.3 our exact result starts to devi- 
ate substantially from the linear order and this deviation 
is particularly strong as the critical value Za = 1/2 is 
approached. The exact A(Za) behavior is also stronger 
than the one found in a recent numerical calculation [20( , 
which also exhibits lattice size dependence. Around the 
critical value Za ^1/2 the function A(Za) has the fol- 
lowing expansion 



A(Za) = 1.12 



\/i - Zcv-0.29 (- - Za 



(11) 

Notice that the induced charge Qind — — | e | ^4 < has 
a screening sign, leading to a decrease of the effective 
impurity charge: Z c q = Z — A(Za). In fact complete 
screening seems possible for Z = 1 and a f=a 1/2, where 
Z e fj = 0. However we will show that the self-consistent 
treatment of the problem (within the Hartree approxi- 
mation) can drastically change this behavior. 

Self- consistent screening. Since the induced charge 
is fully concentrated at the origin, one can easily take 
into account electron-electron interactions in the Hartree 
approximation. The Hartree contribution is expected 
to dominate over the Fock (exchange) contribution for 
TV 3> 1 (N = 4, from the spin and valley degeneracy 
in graphene). To find the effective charge Z e ^ in the 
Hartree approximation, it is sufficient to solve the fol- 
lowing self-consistent equation 

Z c ga = Za — aA(Z c ga). (12) 
If one uses the function A(Za) = A^ = n(Za)/2 calcu- 
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lated in the one-loop approximation, then this equation 
is equivalent to the usual random phase approximation 
(RPA). However, since the exact A(Za) accounts for all 
orders in Za, Eq. (fT2|) is more accurate than the RPA. 

The solution Z c g of Eq. (TT2]) is a function of Z and 
a. The function Z^a versus Za is shown in Fig. [2] for 
different values of the bare charge Z . An impurity with 
charge Z = 1 represents the most important practical 
case. Interestingly, the impurity with Z = 1 remains sub- 
critical for all values of a, i.e. Z e QCt < 1/2. An impurity 
with Z = 2 becomes critical at af =2 = 0.568, and for 
Z = 3 the critical point is af =3 = 0.266. It is clear from 
(fT2|l that the critical for a given bare Z is described 
by the simple formula af = 0.5/(Z — A(0.5)), where 
A(0.5) = 1.12. The fact that A(0.5) > 1 prevents the 
impurity with Z = 1 from becoming critical. Ultimately 
the strongly non-linear variation of A[Za) (Fig. [T]) is re- 
sponsible for the large value of ^4(0.5) and thus the un- 
conventional behavior in the Z = 1 case. 

Even though the supercritical regime is never reached 
for Z = 1, the screening is very substantial, as shown 
in Fig. [^Inset). For example for a — 0.5, Z = 1, we 
find Z e fl/Z = 0.55, and similarly for Z = 2. Thus we 
find quite generally that vacuum polarization screening 
in graphene is very strong, even in the subcritical regime. 
The Z — 1 case is the most relevant experimentally, since 
alkali atoms, such as potassium (K) [10(, typically serve 
as charged scatterers in graphene. For graphene on SiC>2 
substrate with typical dielectric constant e sa 4, the value 
of a is a « 0.9 (using e 2 -> 2e 2 /(l + e), while the RPA 
correction is already taken into account via Eq. (fT2|) ). For 
such large a the vacuum polarization effect is very strong 
and could be important for interpretation of experiments 

To summarize, we have presented an exact solution 
of the Coulomb impurity screening problem in graphene 
in the subcritical regime. The explicit result for the in- 
duced vacuum polarization charge, Eq. is valid to all 
orders in the potential strength Za. The exact solution, 
when combined with the Hartree approximation, leads 
to a self-consistent solution of the impurity problem. An 
impurity with charge Z = 1 is found to be always in the 
subcritical regime, where the induced charge is localized 
at the impurity site. In this regime vacuum polarization 
screening is very strong and weakens substantially the 
impurity potential. However impurities with Z = 2, 3 



and higher can become supercritical. 
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